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Abstract 


Solving the Euler Equations numerically is one of the central issues of the CFD 
community. Accuracy of these solutions not only depend on the scheme used, 
but also strongly depend on the Wall Boundary treatment. We discuss both the 
Strong Fomulation where the zero normal velocity on the Wall is enforced explicitly 
and the Weak Formulation where the zero normal velocity on the Wall is enforced 
through the Flux function. Taking into account the advantages of both these 
Boundary Conditions, Mixed Wall Boundary Conditions have been proposed. We 
discuss here, the Merits and Demerits of the Mixed Wall Boundary Conditions with 
respect to the Mirror Wall Boundary Conditions which is a Weak Formulation. 
Numerical Comparison between the Mixed and Mirror Wall Boundary Conditions 
for the Subsonic and Transonic cases has been made for the Airfoil NACA0012 
with the van Leer and Roe Schemes. The Comparison between the two Boundary 
Conditions has also been made for the Supersonic flow through an Expansion 
Channel. These Comparisons point out that the Mixed Wall Boundary Conditions 
are less dissipative than the Mirror Wall Boundary Conditions. 
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Chapter 1 


Introduction 


1.1 Introduction 


The Engineering fluid dynamic problems of the eighteenth, nineteenth and 
early twentieth centuries almost always involved either the flow of liquids or the 
low speed flow of gases ;for both the cases the assumption of constant density 
is quite valid. However, with the advent of high speed flows, exemplified by 
the de Laval’s convergent - divergent nozzle design and the supersonic sonic 
flight of the Bell Xs-1, the density can no longer be assumed constant through 
out the flow field. Indeed, for such flows the density can sometimes vary by 
orders of magnitude. In this light, we deal here with compressible flows, in 
which the density is not constant. 

The mathematical description of fluid motion in the form of PDE’s by 
Euler not only enriched the science of mathematics but also was responsible 
for the development of a number of useful analytical tools in theoretical fluid 
dynamics as solutions to these PDE’s. A finer mathematical description of the 
fluid, wherein the viscous effects are included was given by Navier and Stokes 
in the nineteenth century. From there on the fluid dynamicists have been 
relentlessly attempting to obtain the closed form solutions to these equations. 
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Computational Fluid Dynamics(CFD) means solving numerically, Euler 
and Navier-Stokes equations either in their original form or approximate form 
on a discretised space or time domain. With the advent of high speed com- 
puters, CFD has gained the credibility as a reliable design tool in the latter 
half of this century. Many numerical schemes have been proposed to solve 
the Euler equations. The schemes of Lax and Wendroff [2], Mac Cormack [3] 
and Jameson [4] are some of the well known schemes, which are purely math- 
ematical in nature, where as Upwind schemes are based on physics of the fluid 
flow. There are two categories of Upwinding schemes such as Finite Vector 
Splitting(FVS) and Finite Difference Splitting(FDS) schemes. The principle 
of FVS is to split flux vector into a positive and negative parts such that the 
positive flux Jacobian has all the positive eigenvalues and the negative flux 
Jacobian has all the negative eigenvalues. The FDS schemes are based on the 
classical initial value problem called Riemann problem. The accuracy of the 
solutions thus obtained not only depend on the accuracy of the scheme, but 
also strongly depend on the Wall boundary treatment. The solution accu- 
racy and stability very strongly depend on the boundary procedures adopted 
in enforcing the wall boundary conditions. So, proper treatment of the wall 
boundary conditions is essential for obtaining reliable results. Therefore, we 
have tried to look into this aspect here. 


1.2 Classification 

Accurate satisfaction of the boundary conditions is one of the central issues in 
any CFD problem. A solution acquires a special flavor because of the bound- 
ary conditions of the problem. One way of classifying the boundary conditions 
is (a). Physical boundary conditions such as flow tangency, no-slip, pressure on 
the exit plane etc. (b). Fictitious boundary conditions which are usually re- 
quired at the farfield boundary of a finite sized computational domain. In the 
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language of Moretti such farfield boundary conditions are mathematical model 
for the rest of the universe, van Leer calls (a) type Boundary Conditions as real 
Boundary Conditions and (b) type Boundary Conditions as artificial Bound- 
ary Conditions. There is one more way of classifying Boundary Conditions, 
(a) Dirichlet Boundary Condition, (b) Neumann Boundary Condition. If the 
Boundary Conditions are directly specified on the unknowns u and v, then 
such Boundary Conditions are Dirichlet Boundary Conditions, where as if the 
Boundary Conditions are specified on the derivatives of the unknowns, such 
type are called Neumann type of Boundary Conditions. For example, no-slip 
condition for viscous flows is a Dirichlet Boundary Condition while, the flow 
tangency condition in terms of the full potential is a Neumann type Boundary 
Condition. 


1.3 Literature Survey 


Two basically different approaches can be found in the literature to introduce 
Wall Boundary Conditions for the Euler equations. In one approach, the 
unknown along the boundary are determined from inside the domain by some 
space extrapolation technique under the constraint that the unknowns must 
satisfy the Boundary Conditions. Various methods have been proposed. The 
use of normal momentum equation for pressure extrapolation at a solid wall 
belongs to this approach. 

A second approach follows more closely the physics of the hyperbolic 
equations by discretising the compatibility equations for the outgoing charac- 
teristic variables, while replacing the ingoing characteristic information by the 
physical Boundary Conditions. Consistent Boundary Conditions for the cell 
centered upwind finite volume Euler solvers proposed by Deconinck and Struys 
[5] belongs to this second approach. In this approach, from the consistency 
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condition for a numerical flux function, the unknowns along the boundary are 
determined such that the numerical flux function based on the boundary and 
adjacent interior variables leads to a flux which satisfies the Boundary Condi- 
tion. This boundary treatment reduces to characteristic boundary conditions 
for the case of a linear hyperbolic systems of equations. 

Direction of information propagation at the boundaries must be properly 
taken care. As Butler [7], points that use of a special scheme at the bound- 
ary, violates the global conservation, causes inconsistency, and many lead to 
artificial boundary layer of numerical error. Instead he proposes to evaluate 
boundary points by relating them to interior values by the theory of character- 
istics. This approach has been used for treating conditions at both shock-wave 
and solid boundary [8]. An alternative to the theory of characteristics at a 
solid body is to use the knowledge that a streamline in two dimensions and 
a stream surface in a three dimensions wet the surface of a body in inviscid 
flow. This information, when coupled with the solid wall boundary condi- 
tion and incorporated into the momentum equation, which leads to relation 
that connects boundary points to those in the interior. For sometime various 
workers used this so-called normal momentum equation but usually for only 
simple body shapes and coordinate systems. Mac Cormack [9], presented a 
general expression valid for any geometry, but one that requires labourious 
interpolation between grid points. Later on Rizzy [10], gave an equally general 
formulation that avoided the need for interpolation. Followed by that, Rizzy 
[11], presented the derivation for his auxiliary equation and demonstrated how 
it is computationally more convenient to apply because it requires no inter- 
polation of computed flow properties, but rather uses data only at the grid 
points, a particular advantage in three dimensions. 

One of the common procedures found in the literature is the Mirror con- 
dition in which the approximate Riemann solver used for determining fluxes in 
the interior edges is used on the wall boundary also. For this, a fictitious state 
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is determined in such a way that the zero normal mass flux across the wall is 
guaranteed. So, based on the information obtained from interior and fictitious 
states, the fluxes on the wall are calculated in a way similar to that used in 
the interior. Its main advantage is that it takes into account, the direction of 
propagation of information.lt also converges fast. Its main disadvantage is the 
numerical dissipation that it produces. Moretti [12] is a strong critic of this 
method. While criticising this method he says that many authors use what 
they call a reflection principle which is not a principle at all and should rather 
be called a reflection technique. He also points out that there is no precise def- 
inition for such a technique in the literature. All the Physical parameters are 
specularly reflected on the rigid wall except the normal velocity which is as- 
sumed as antisymmetric. These assumptions force the normal derivatives of 
these variables to vanish at the wall. Here, the assumption = 0 doesnot 
take into account the curvature effects. Hence, he argues that even some of 
the assumptions made in this procedure are physically wrong. 

Various approaches can be found in the literature to introduce wall 
Boundary Conditions for the Euler equations. But, all these schemes one 
way or the other way boils down to the enforcement of zero normal velocity 
on the solid wall boundary. Based on the method of enforcement of zero nor- 
mal velocity on the solid wall, these formulations can classified into strong 
and weak formulations. In the strong formulation, the zero normal velocity on 
the wall is enforced explicitly, where as in weak formulation, the zero normal 
velocity is enforced through the flux function. Taking into advantage of both 
the approaches Mixed boundary procedures have been developed [1], 

Improper treatment of wall boundary conditions allow numerical dissi- 
pation to creep into the solution. So, the Wall Boundary Conditions play an 
important role in obtaining the proper results. The Strong Formulation which 
we discuss later, violates the global conservation whereas, the Weak Formula- 
tion such as Mirror Wall Boundary Condition donot take into account the cur- 
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vature effects and are less dissipative. The Mixed Boundary Condition which 
has the advantages of both these Schemes is expected to be less dissipative. 
Since, the numerical dissipation can spoil any solution, it becomes important 
to keep this in check. It is in this context that we are looking at the Mirror and 
the Mixed Wall Boundary Conditions which are discussed in the subsequent 
chapters. In our work, we have applied the above stated Mixed boundary 
procedure to cell centered finite volume schemes on structured meshes. Here, 
we have tried to examine the problem of numerical dissipation in Mixed and 
Mirror wall boundary conditions. 



Chapter 2 


Formulation and Flux formulae 


2.1 2-D Euler equations 


The 2-D Euler equations of Gas Dynamics can be written in the differential 
form of conservation law as, 


dU dF dG_ 
dt dx dy 


( 2 . 1 ) 


where, U is the vector of conserved variables and F and G are flux vectors in 
the x and y directions respectively. They are defined by the relations: 



“ p 


pu 


’ pv 

u = 

pu 

;F = 

pu 2 + p 

; G — 

puv 


pv 

puv 


pv z +p 


e 


. (e + p) u . 


_ (e + p) v m 


where, p is the mass density, u is the x-component of the velocity vector, 
v is the y-component of velocity vector, p is the pressure and e is the total 
energy per unit volume given by, 


p P( u f\ + u \) 

7 - 1 + 2 


e = 


( 2 . 3 ) 
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2.2 Finite Volume Schemes 

Using the Gauss theorem the state update formula for finite volume is written 
as 

{ift—kLww* (2 - 4) 

where, the vector of state variables U = (p, pu, E) T represent the mass density, 
momentum components and total energy per unit volume. The subscript i 
represents the i th finite volume The Euler flux normal to the finite volume 
interface dQ l is given by 

F± = F.h = (pu_ L , pu\u + pn, (E + p)u±) T (2.5) 

where, h is the outward unit normal, ~ u.h and My = uAn are the normal 
and tangential velocities respectively. The pressure p is given by 

P = (7 - 1) ~ 1 + «|)) (2.6) 

where, 7 = 1.4 is the specific heat ratio. 

It is possible to build finite volume schemes either in cell vertex or cell 
center framework. Typical cell vertex and cell center finite volumes on a quadri- 
lateral domain are shown in fig. 2.1: 


(a) Cell Vertex. _fra.meioorlc (b) Cell center Jrameiucrlc 
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Figure 2.1: Control volume 


In cell vertex finite volume framework, the state variables are updated at 
the nodes whereas, in cell center finite volume framework, the state variables 





2.2 Finite Volume Schemes 


9 


are updated at the centroids of the cells. We shall refer to these points as 
reference points. 

In the case of cell center finite volume schemes, as the reference points 
lie at the center of the cell, they donot lie on the boundaries. Whereas, in the 
case of the cell vertex finite volume schemes, as the reference points are the 
nodes of the cells, they lie on the boundaries. So, there will be difference in the 
way the boundary conditions axe implemented in these two schemes. In cell 
vertex finite volume schemes, the boundary conditions axe enforced explicitly 
on the reference points whereas, in the cell center finite volume schemes, the 
boundary conditions are enforced through the flux function. As we have used 
the cell center finite volume schemes, we have applied the latter method. 


2.2.1 Higher Order Finite Volume Schemes 


For steady state computations, the accurate calculation of fluxes on the cell 
interface using higher order spatial accurate schemes is very important. The 
time integration can be done by using either implicit schemes or multistage 
schemes. For achieving high order spatial accuracy on structured meshes, we 
are made using the MUSCL (monotonic upstream centered schemes for con- 
servation laws.) ( Anderson, et.al [13] ) extrapolation. The conserved variables 
axe extrapolated on to the cell interfaces for subsonic and supersonic cases 
and for hypersonic cases, the primitive variables (as done by Sreekanth and 
Reddy [1992a] are extrapolated. The one parameter family of schemes called 
k - schemes are are used in MUSCL extrapolation are given below : 


W L = W k + &- ((1 - k) (W k - W k ^) + (! + *) (W k+1 - W k )) ' 
W R = W k -& ((1 + k) (W k+1 - W k ) + (1 - k) (W k+2 - W k+1 )) , 


(2.7) 


where W represents either the conserved variable U or the primitive variable 
V, as the case may be. For <f> = 0 we obtain a first order scheme. For <j> = 1 
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where W represents either the conserved variable U or the primitive variable 
V, as the case may be. For <j> — 0 we obtain a first order scheme. For <£>= l 
we have the following schemes for different values of k: k = — 1 (second order 
- Upwind), k = 1 (second order - Central ), k = 0 (Fromm Scheme), k — | 
(QUICK), k = | (Third Order - Upwind Biased ). The spurious oscillations 
that appear in the solution because of use of high resolution schemes can be 
suppressed using van Albada limiter ( Anderson, et.al [13] ): 


2A + W k A„W k + e 
A + W k 2 + A_W k 2 + e 


(2 8 ) 


where, 

A+W* = W k+l -W k \ 

A -W k = W k -W k - 1 j 

Here e is a small number added to avoid fault of divide fault during the com- 
putation. In all the higher order computations on structured mesh presented 
in this thesis k is choosen to be |. 


2.2.2 Cell calculations 


The space discretised finite volume form of the Euler equations is given by, 

= ( 2 . 10 ) 

where the subscript k denotes the cell number, m is the edge index (or interface 
index) of the k th cell and the subscript n denotes the time step. F m and As m 
are respectively the normal flux and the area of the m th face. Further, V k is 
the volume of the k th cell. The vector of conserved variable U and the flux 
vector F m are given by, 


p 


' pu L 

pu 

■ FT = 

puux + pn x 

pv 

i x j. 

pvux + pn y 

e 


. (e + p)ui . 


( 2 . 11 ) 
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Where, n x and n y represent the direction cosines of the unit normal h to the 
interface and are given by, 


Ay 

n * = A = 


A.x 

As 


( 2 . 12 ) 


It is very essential to know about the cell and its related data for the under- 



X 


Figure 2.2: A typical two dimensional finite volume cell. 


standing of any numerical code. So we are presenting a typical finite volume 
cell and the way fluxes are calculated on each edge of the cell as shown in the 
fig 2.2. 

2.2.3 Method of Least Squares 

The gradient terms in this work are obtained using the method of least squares. 
Here the gradient terms are obtained by minimising sum of the squares of the 
errors [14], defined by 

E, = £ (A*. - (V«i • Arj) 2 (2.13) 

J=l 

where, subscript j represents the neighbouring grid points, N is the number of 
support points, 4> t is the value obtained by the state update formula, = 
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Figure 2.3 A typical 2-D Finite volume Cell and its neighbours. 


4>j — (j) t and Afj — r \ — r, We have considered all the neighbouring points ie 
N — 8 as shown in the fig2.3. 


For a linear reconstruction procedure, it is possible to obtain a closed 
form expression for the gradient terms by minimising E t which is defined above 
(eqn.no 2.13) and they are given by, 

|| Ay || 2 (Ax, A(J>) - ( Ax , Ay) (Ay, A<j>) 


4*i, x — 

4>i,y ~ 


II Ax || 2 || Ay || 2 - (Ax, Ay) 2 

Ax || 2 (Ay, Afi) - (Ax, Ay) (Ax, A^) 
II Ax || 2 || Ay || 2 - (Ax, Ay) 2 


(2.14) 

(2.15) 


where, 


II Ax || 2 = Yj Ax j> (A^, A<f>) = E^=i AxjAcfjj 

j = l 

II Ay H 2 = E (Ay,A« = E",A ft A^ 

(2.16) 


The above expressions (3.. 14- ) and (US’) for the determination of the gra- 
dients of the variables are used in satisfying gradient boundary conditions 
on both structured and unstructured meshes. To calculate the Pressure at the 
Wall satisfying the normal momentum equation, the gradients of the conserved 
variables are required, these gradients are calculated using the above stated 
Least Square Formulae. 
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2.3 Flux Calculation 

2.3.1 Hyperbolicity of Euler equations 


One of the important features of the Hyperbolic conservation laws is that they 
have preferential direction of information propagation. This property is called 
as hyperbolicity property. This can be best understood by casting the Euler 
equations in the form of scalar convective equations. The lD-Euler equations 
of Gas Dynamics can be written in the differential form of conservation laws 


as 


dU dF_ 
dt dx 


0 (2.17) 

where, U is the vector of conserved variables and F is the flux vector. They 
are defined by the relation 


U = 


~u x 1 


* p 


pu 

u 2 

= 

pu 

•F = 

pu 2 +p 

U 3 


e 


(e + p)u 


(2.18) 


Where, p is the mass density, u is the fluid velocity, p is the pressure and e is 
the total energy per unit volume given by, 

,2 


e = 


p pu 

(WT) + 1" 


(2.19) 


Eventhough F (2.18) is expressed as a function of primitive variables, p,u,p it 
can also be expressed in terms of the components of the vector of conserved 
variables U. We then obtain 

V 2 


F = 


( 7 - m + 


I v. 


( 7 - 1 ) V{ 

2 m 


( 2 . 20 ) 


Using the fact that F = F(U) the equation(2.20) can be cast in the form 


of 
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BU * dU n 

~m +A a^~° 

where the flux Jacobian matrix A = §§ is given by 


(2.21) 


A = 




(3-7)11 (7-1) 


. g2u _l iirJl u 3 ... °L 

(7-1) + 2 u (7-1) 




-yu 


( 2 . 22 ) 


Here a is the sonic speed given by a = J (yp) / p. The matrix A has all real eigen 
values and a complete set of linearly independent eigen vectors. Therefore the 
Euler equations are non-linear hyperbolic PDEs for the unknown U and must 
be solved with suitable boundary conditions. To cast Euler equations in the 
form of scalar convection equations, the flux Jacobian A is cast in the Jordan 
canonical form as follows: 


A = RAL = RAR- 1 


( 2 . 23 ) 


where, R and L are the matrices of right and left eigen vectors of A respectively 
and A is a diagonal matrix, the entries along the diagonal being the eigen values 
of A. The matrices R,A and L are given by: 


1 1 1 


u — 

a 

u 

u 

+ a 

it 2 1 a 2 

2 -««+( 7-1) 

V? 

2 

\ + ( 7 - 1 ) 

JL _L tzD. u 2 - 

2 a + 4 a 2 U 
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2a “ 

_ IzzAL b- 1 ) 

2 X U 2a 2 

1 _ 

il 
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tl zll 

a2 

_jl 4. 2 

2 a + 4 a 2 U 
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£cii„ 

2a 2 U 

hzH 

2a 2 


u — a 

0 
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A = 

0 

u 

0 



0 

0 

u + a 



( 2 . 24 ) 


( 2 . 25 ) 


( 2 . 26 ) 
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Therefore, substituting (2.20) into (2.18) gives 


au „, r au „ 
aF + SA£ aF = 0 

(2.27) 

Multiplying by L on both the sides, we get 


T dU kT dU 

£ ir +A£ & =0 

(2.28) 

Letting LdU = dip, we get 

dt +A dx ° 

(2.29) 


It should be noted that the convection equations defined above, are coupled 
through the dependence of A, and ip t on the components conserved variable 
vector U. Under the isentropic conditions, the t/>,’s reduces to the well known 
riemann invarients. The hyperbolicity of Euler equations become more explicit 
and clearer from thier convective form of equations (eqn.no. 2. 29). The quan- 
tities ip t get convected with the speeds of the eigen values of the flux Jacobian 
A. Numerical schemes which follow the physics of the flow should take care 
of the propagation of information. The recent Upwind schemes such as Flux 
Vector Splitting and Flux difference Schemes take care of this hyperbolicity 
property of the Euler equations. 

2.3.2 Flux Vector Splitting Schemes 

Flux vector splitting is a technique for splitting the flux vector into upwind 
components . Here, the flux vector F is split into a positive part F + and a 
negative part F ~ , such that the positive flux Jacobian has all the positive 
eigenvalues and the negative flux Jacobian has all the negative eigen values as 
follows : 

F(U) = F + (U) + F-(U) (2.30) 

where , F + and F~ are are constructed such that 

dF + 

— ([/) = A*(U) 


(2.31) 
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Figure 2.4: Flux vector splitting 


-wW-a-V) 

The above equations imply that the eigenvalues of A + are nonnegative (> 0) 
and those of A ~ are nonpositive(< 0). The numerical flux formula at an 
interface is given by 

F(U l ,Ur) = F + (U l ) + F~(U r ) 


By applying FVS, the Euler equations become 


dU dF+_ dF^ 
dt dx dx 


= 0 


(2.32) 


In the above equation the positive flux is approximated using with the back- 
ward difference operator and the negative flux is approximated with the for- 
ward difference operator. There is no unique way of splitting the fluxes. We 
discuss one such FVS schemes due to van Leer [15] which we have used in our 
work in the subsequent sections. 


2.3.2. 1 van Leer scheme 

van Leer has proposed a method of splitting the flux vector [15]. The mass 
momentum and energy fluxes were split such that the forward and backward 
fluxes transitioned smoothly at eigenvalue sign changes ie near sonic and stag- 
nation points. The flux F(W) with W = (p, a, M) T is split into a forward flux 
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F + (W) and a backward flux F (W) as follows: 

FJW) = Fl(Wn) + FZ(W l ) 
with the local algeobraic mach number — uj_/a , 

1 


n = + ‘) 2 




~^(M ± - 1) ! 


a £ ld -yj.+2«l _ nyU{l 

n^y-l^+ia) + nxU|| 
[{~i~l)u L +2a] 2 +u\ 


»- ]_ WbU|| 

”y[(7-l)^x-2a] _j_ n 


U|| 


[(•Y-l)u x -2a] 2 +u2 


(2 33) 


(2.34) 


(2.35) 


2.3.3 Flux Difference Splitting(FDS) 

They fall into the category of solving the approximate Riemann solvers. In 
FDS schemes, the flux difference is split into upwind, components. Here, the 
flux difference (A F) is split into a positive part and a negative part such that 
the positive part contains all the positive eigenvalues and the negative part 
contains all the negative eigenvalues . 

A F = (A F)~ + (AF) + . (2.36) 

As in the FVS schemes, there is no unique way splitting the flux differ- 
ence.Different schemes have been proposed based on the method of splitting . 
we discuss one such method due to Roe[16] which is used in this work. 
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2. 3. 3.1 Roe scheme 

Roe has proposed a FDS scheme which is named after him [16]. According to 
his scheme, 

dF = AdU (2.37) 

The finite difference analogue of the above differential relation is , 

A F = AAU (2.38) 

where, 

A F = F r - F l (2.39) 

A U = U R - U L (2.40) 

In the above equation subscripts R and L represent the right and left states 
respectively. Roe[16] has given the methodology for the construction of A 
called the Roe averaged matrix. Expressing A in the Jordan canonical form, 

A F = RALAU (2.41) 

By introducing A# = LAU , 

A F = RA AV (2 42) 

or, 

A F= J2 A.r.Atf, (2.43) 

t=l,3 

Now the term A F is split into two parts (A F)~ and (A F) + , which are the 
fluctuations corresponding to the backward moving waves and the forward 
moving waves respectively. This is done as follows: 

(A F)~ = £ \ t r t A% (2-44) 

A ,<0 

and, 

(AF)+ = Y. 

>>.>0 


(2.45) 
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Therefore it is possible to calculate the flux at the interface Fj using one of 
the formulas given below: 


F T = F L + (A F)~ (2.46) 

Fr = F R - (A F) + (2.47) 

Generally the the interfacial flux is expressed in terms of the symmetric formula 
obtained by taking average of the fluxes obtained using the above two formulas. 
That is, 

Fi = Uf l +Fr )- 1 ~y: a i a ‘ i ( 2 - 48 ) 

1 1 t=l,3 


The flux function of Roe [16] reads, 

F±(W r ,W l ) = ^(F ± (W r ) + F x (W l ))~ i I A(W R ,W L ,n) I (W R -W L ) (2.49) 
where, A{W R , WL,h) is the Jacobian matrix defined using Roe’s average: 


y/PRPR + y/PLPL 

(2.50) 

y/pR + yfpL 

y/PRU R + yfpLUL 

(2.51) 

y/pR + y/pL 

y/pR.H R + y/plHt 

(2.52) 

yfpR + \[PL 


2.4 Summary 

In this Chapter, we have discussed the basic ideas of Finite Volume Schemes, 
The Least Square Method, which has been discussed is useful in obtaining 
the gradient terms. Two important methods of obtaining the fluxes on the 
Cell interfaces such as Flux Vector Splitting Scheme and the Flux Difference 
Scheme have been discussed. The van Leer Scheme which is a FVS Scheme 
or the Roe Scheme which is a FDS Scheme is applied in conjunction with the 
Mirror or Mixed Boundary Condition in the various test cases that are carried 
out. These Boundary Conditions axe discussed in the subsequent Chapters. 



Chapter 3 


Wall Boundary Conditions 


A lot of work has already been done on finding of fluxes for the interior edges. 
But not much attention was given to the Wall Boundary Conditions. For the 
Computation of many Fluid Dynamic problems more difficulty is encountered 
in satisfying the Boundary Conditions than in balancing the Differential Equa- 
tions at the interior points of the flow field. This is so because on the Boundary 
not all the flow Variables are specified by the Boundary Conditions, and there 
remain more unknowns than the interior Equations. The problem is fully de- 
fined only when the proper Boundary Conditions are specified. The Boundary 
Conditions are such an important part of the Problem that the pattern of 
the two flow fields can be completely different from one another because of 
some differences in the flow boundaries, despite the fact that both Flows obey 
the same system of Partial Differential Equations. Flow essentially develops 
its characteristics because of what happens on the Boundaries. Thus Wall 
Boundary Conditions becomes important. 

From the theory of characteristics, at any boundary, the variables corre- 
sponding to the outgoing and incoming information must be identified. Only 
the variables corresponding to the information transported from the boundary 
towards the interior ^incoming') can be freely imposed as physical boundary con- 
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ditions. All other variables which correspond to the outgoing information must 
be evaluated from the interior information which are called numerical bound- 
ary conditions. As for as the solid wall is considered, only one characteristic 
enters the flow domain from outside and the rest of the characteristics leave 
the flow domain. Refer fig 3.1. Therefore it is necessary to provide only one in- 
formation on the wall based on the physical considerations as natural boundary 
conditions. This condition is expressed by the vanishing normal velocity. 

u x = 0 (3.1) 

As a result all the convective components through the solid wall will vanish in 



Figure 3.1: Directions of information flow at the Wall boundary 


the computation of the flux terms and the flux vector reduces to the following 



3.1 Strong formulation 


22 


expression in a 2-D flow. 


0 


Fxb 


pri x 

pn y 


0 


(3.2) 


Hence, only the pressure component remains at the wall. So, in the Finite 
Volume schemes, only Pressure is required for calculating Fluxes on the Wall 
Boundary Conditions. 


As we stated earlier, there are two ways of implimenting the Wall bound- 
ary conditions. One of them is the Strong formulation while the other is the 
weak formulation. In the following section, we present the strong formulation 
in detail which is followed by the weak formulation. 


3.1 Strong formulation 


In the strong formulation, the zero normal velocity on the wall is enforced ex- 
plicitly, ie state variables on the wall boundary are defined using the physical 
constraints which are explained in this section later, where as the state vari- 
ables at the interior nodes are updated using the Finite Volume state update 
formula(eq.no 2.4). 

As we stated earlier, at the wall boundary, only one information has to 
be provided explicitly at the wall. Naturally this condition will be 

ux = 0 (3.3) 


The second physical constraint corresponds to the normal momentum 
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equation for the rigid stationary walls satisfying = 0 and is given by 


Tn “ - p R 

where p is the mass density, ity is the tangential velocity and R is the radius 
of curvature. But the use of the radius of curvature in the above equation 
can be avoided simply by using the differential relation obtained from the non 
conservative form of Euler equations. 

dp 

— = -pn.((u.v)u) (3.5) 

It is also possible to obtain a similar relation for the normal gradient of pressure 
on the Wall using the conservative form of Euler equations. One more way of 
determining the pressure on the wall satisfying the normal momentum equation 
is to use the conservative integral relation [20], 


F\ ds 


Jdi 


Fids 


Here, VLb and fl/ represent the boundary and the interior edges of any wall 
boundary cell respectively. Assuming that the pressure is constant along the 
boundary edges of each cell of the wall boundary cells, we have the following 
relation for the pressure on the boundary, 


Pb = 


[ F X2 ds + n f F ±3 ds 
Jd n 7 Jan T J 


where, Asg is the length of the boundary edge. In our work we have used (eqn 

no.3.7) for satisfying the normal momentum equation on the wall. The third 

physical constraint used, is the Crocco’s relation [19]. Crocco’s relation in a 

direction normal to the wall is given by 

dS 7 - 1 f oN 

51 = -—” 11 “ < 3 ' 8 ) 

where , S = p/p 7 is an entropy like variable. Here isentropy of the flow is 

assumed. In the case of inviscid subsonic flows without discontinuities, this 

relation becomes, 

^ = 0 (3.9) 
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The computations in this work are based on the (eqn.no. 3. 9). Enforcing the 
normal entropy gradient on the wall boundary has been found less dissipative 
for subsonic and transonic flows. The computations in this work, constraint 
(eqn.no.3.9) is used in preference to (eqn.no. 3.8). Even if the flow is irrota- 
tional(as in the subsonic test case), the vorticities computed on the wall are 
considerably high [1] because of the numerical errors. This problem is found 
to be even more accute in the stagnation regions. Therefore the use of such 
erroneous wall vorticities in the satisfaction of the entropy condition on the 
wall makes the boundary condition more dissipative. The fourth constraint 
used is related to the total enthalpy, ie Total enthalpy remains constant for 
the inviscid flows without source terms if the flow at infinity is isenthalpic. ie 


dH 

dn 


(3.10) 


where, H = ( E+p)/p is the total enthalpy. Thus the above equation can alsi 
be used as a constraint for obtaining the boundary conditions. 


The values of any variable $ (S,H or p as the case may be) on the wall 
can be obtained by enforcing the conditions on the normal gradients using the 
Method of Least Squares ie (eqn. no. 2.13). 

The main advantage of the strong formulation is that it is less dissipa- 
tive. It also takes into account, the curvature effects of the wall boundary. The 
curvature effects axe taken care by the equation number 3.4. The main disad- 
vantage of the strong formulation is that, the scheme used for the boundary 
points is different from the one that was used for. the interior points. There- 
fore it violates global law of conservation. The other drawback of the strong 
formulation is that, its implementation is direct in the case of cell vertex finite 
volume framework where as in case of the cell center finite volume framework, 
it is not so. This is because, the reference point will be present on the bound- 
ary surface in the case of cell vertex framework where as, it is not so in the 
case of the cell center framework. Also in the case of the cell center framework, 
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the normal velocity is not zero at the reference points of the boundary cells. 
Hence, (eqn.no.3.1) is not valid in this case. For structured meshes this prob- 
lem can be overcome by using the half cells at the wall boundary. Therefore it 
is necessary to design weak boundary conditions which take into account the 
curvature effects and also also which are less dissipative. 


3.2 Weak formulation 

In the weak formulation, the natural boundary condition = 0 is enforced 
through flux function. Therefore the weak boundary conditions reduces to 
the determination of fluxes at the cell interfaces coinciding with the bound- 
ary through the physical constraints. The reference points of the boundary 
are updated in a way similar to the interior cells. So contrary to the strong 
formulation, the implementation of the weak boundary conditions is same in 
both the cell center and cell vertex frameworks. Here, we are presenting Mirror 
boundary condition which we have used in our work which is later followed by 
the Mixed Boundary Condition in next section. 


3.2.1 Mirror condition 


In this boundary condition, the variables in the fictitious cells are defined 
such a way as to ensure the vanishing mass Flux at the wall. Once, the 
fictitious state are defined, the fluxes on the boundary edges are calculated 
based on the information obtained from the computational domain(interior 
state) on one side and the fictitious cells on the other. So, in this boundary 
condition the approximate Riemann solver used for determining the fluxes on 
the interior edges is used on the wall boundary also. There are various ways of 
defining the fictitious state. In our work we are using the reflection principle 
which guarantees the zero normal flux across the wall boundary. As the name 
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Figure 3.2: Mirror boundary condition 


indicates the pressure, density and the tangential velocity of the interior cells 
are reflected on to the fictitious cells. To make the mass flux on the wall to be 
zero, the normal velocity in the fictitious cell is equal and opposite to that in 
the interior cell. Therefore the mirror condition is stated as follows: 

U_Lf = —u±j 

u \\f = u ¥ 

( 3 . 11 ) 

Pf - Pi 
Pf = Pi 

where, the subscripts I and f refer to the interior and fictitious states respec- 
tively. 

The important advantage of the Mirror condition is that it permits the 
use of the interior scheme and thus the direction of information propagation 
is maintained at the boundary also. The main drawback of this boundary 
condition is that it is dissipative and it doesnot take into account the curvature 
effects. 
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3.3 Mixed Boundary Condition 

It is clear from the above discussion that both the strong formulations and the 
weak formulations have their own merits and demerits. In order to harness the 
advantages of both the formulations, we are using the Mixed boundary condi- 
tion. One of the ways incorporating the features of the Strong Formulation in 
the Weak Formulation is by making use of the constraints in the Equations 
numbered 3.3, 3.7, 3.8, 3.10 directly to obtain the wall fluxes. But such a 
procedure boils down to determining simply Pressure on the wall using the 
normal momentum equation. This is found to be unstable and leads to the 
development of artificial boundary layer on the wall. Therefore it becomes 
important to look into other ways of implementation and Mixed Boundary 
Condition is a result of such an outcome [1]. 

In this Mixed boundary treatment, we define the fictitious cells by using 
the strong formulation, that is by using the physical constraints stated in the 
strong formulation. The interior state is determined in a way similar to the 
Mirror Condition. Now both the states across the boundary edges are defined 
and so the approximate Riemann solver can be used to determine the wall 
fluxes. 

It can be noted here that this boundary treatment is also not conserva- 
tive, that is it allows the transpiration of the mass flux across the rigid wall 
boundary. But, the justification for the use of this treatment comes from the 
fact that the computed fluid velocities normal to the boundary wall at the 
steady state are very small. Numerical errors involved are also very small. 
Numerical results also confirm this fact. 
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3.4 Summary 

In this Chapter, we have discussed the basic ideas of the different types of for- 
mulations used for the Wall Boundary Conditions. In the Strong formulation 
zero normal velocity on the Wall Boundary is enforced explicitly whereas, in 
the Weak Formulations, the same is done through the flux function. As an 
example for the Weak Formulation, we have also discussed about the Mirror 
Boundary Condition. Both these Schemes have their own advantages and the 
disadvantages. The Mixed Boundary Condition which takes into account the 
advantages of both the Strong and Weak Formulations, is expected to produce 
better results. Many numerical experiments have been performed to verify 
this fact. The results obtained with the Mirror Wall Boundary Condition have 
been compared with those obtained with the Mixed Wall Boundary Condi- 
tion which is a Weak Formulation. These results have been discussed in next 
Chapter. 



Chapter 4 


Results and Discussions 


The boundary conditions discussed in the earlier chapters have been tested for 
various test cases such as, 

1. Subsonic flow past NACA 0012 

2. Transonic flow past NACA 0012 

3. Supersonic flow through an expansion channel 


In the first two cases, van Leer, Roe and the Combination of both the 
schemes .i.e., Roe scheme for the interior cell and the van Leer scheme for 
the boundary cells have been used. The Mixed Boundary Condition used in 
conjunction with the FDS doesnot seem to improve the quality of the results 
significantly [1]. Therefore a combination of FDS scheme for the interior edges 
and the FVS scheme for the boundary edges has been carried out in this work. 

The details regarding the grids that are used for the various computations 
axe presented in the table 4.1. All the grids are symmetric regular quadrilateral 
grids. These grids are shown in fig.4.1. In order to extract the information at 
the walls properly, the number of grids at the wall have been increased. In the 
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results presented, the wall entropy deviation is defined as, 



(4.1) 


All the computations were made using either van Leer or Roe’s Scheme. 


Table 4.1 
Grid details 


Name 

configuration 

Nodes 

cells 

Edges 

Gridl 

Airfoil NACA0012 

2112 

2048 

4160 

Grid2 

Expansion channel 

2821 

2700 

5520 


4.1 Subsonic flow past NACA 0012 

This is one of the standard GAMM [16] test cases for evaluating different 
schemes for the subsonic flows. The free stream Mach no is 0.63 and the angle 
of attack is 2.0. 

Mach contours obtained by using van Leer Scheme are shown in fig4.2 
The wall entropy generated for the same case are presented in the fig.4.3. The 
entropy layer generated by the wail boundary condition becomes apparent 
when we look at the Mach contours. It is also clear from the wall entropy 
plots. The wall entropy generated in the case of Mirror Boundary Condition 
is in the range of 0.01 to 0.04 whereas that in the case of Mixed Boundary 
Condition, it is in the range of 0.005 to 0.02. So, clearly the entropy generated 
bv the Mirror Boundary Condition is higher than that generated by the Mixed 
Boundary Condition. 

The Mixed Boundary Condition used in conjunction with the FDS Schemes 
such as Roe, doesnot seem to improve the results significantly over the Mirror 
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Boundary Condition [1]. In order to extend the applicability Mixed boundary 
condition for the FDS schemes, we have tried using the FDS schemes for flux 
calculation of the interior edges and FVS schemes for flux calculation of the 
boundary edges. The Mach contours are shown in fig.4.4. The entropy gen- 
eration at the wall is plotted in the fig. 4. 5. As seen in the earlier case, the 
Mach contours obtained using the Mixed Boundary Condition are superior to 
those obtained using the Mirror Boundary Condition. Again, as seen earlier, 
the Mach contours in this case also clearly shows the superiority of the Mixed 
Boundary Condition over the Mirror Boundary Condition. This is again con- 
firmed by the wall entropy changes which show that the wall entropy in the 
case of Mirror Boundary Condition are much higher than that in the case of 
Mixed Boundary Condition. 

The above results show that the entropy generated on the wall due to 
Boundary Condition in the stagnation region was high irrespective of the type 
of wall Boundary Condition. The entropy generation due to mirror Boundary 
Condition was high as compared to that of the Mixed Boundary Condition. 
In the stagnation region the gradients with the flow field hifh. As a result, the 
numerical dissipation will also be more in this region. So, the Vortices gen- 
erated i this region are also high. These Vortices are convected downstream 
towards the leading edge. To reduce the Vorticity generation in the stagnation 
region, either a fine grid or a numerically less dissipative wall boundary condi- 
tion should be used. By following the latter method, we have used the Mixed 
boundary condition which is less dissipative in the stagnation region and the 
Mirror boundary condition downstream. 

The mach contours and the wall entropy obtained are shown in fig.4.6. 
The Mach contours shows a clear improvement over the Mirror conditions, 
though may not be better than the Mixed Boundary Condition of the van 
Leer scheme. The wall entropy plot confirms this fact much more clearly. 
But, the experiment produces oscillations in the wall entropy generated at the 
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stagnation region. This is not desirous. This is mainly because of the sudden 
switch from the Mixed wall Boundary Condition to the Mirror wall boundary 
condition which are not smooth. 


4.2 Transonic flow past NACA 0012 

In this case, the Mach number is 0.8 and the angle of attack is 1.25. This is 
one of the standard cases of AGARD [17]. This is a very interesting test case, 
because this has a shock in the upper surface and a relatively weaker shock, 
(closer to the leading edge than the upper shock) on the lower surface of the 
airfoil. It is generally possible to capture the upper surface shock accurately, 
but capturing the lower surface shock is difficult(AGARD[17]). The results of 
the transonic flow cases follow the lines of the subsonic flow which were shown 
in the earlier sections. 

The Mach Contours obtained by using van Leer Scheme are presented 
in fig. 4. 7. and the wall entropy plots are presented in fig. 4. 8. This fact is 
also confirmed when we look at the wall entropy plots which show that the 
entropy generated at the wall (which is purely numerical) due to the Mirror 
wall boundary condition is high as compared to that due to the Mixed wall 
boundary condition. 

Similar to the subsonic flows, to get the improved results for the Roe 
scheme , computations were done using the van Leer Scheme at the wall and 
Roe scheme in the interior cells. The results as expected, showed improvement. 
The Mach contours, and the wall entropy plots are shown in fig-4.9 and fig.4.10 
respectively. Both the plots confirm the superiority of Mixed wall boundary 
condition over the Mirror wall boundary condition. 
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4.3 Supersonic Flow Through an Expansion 
Channel 

In this case the free stream Mach number is 2.0. This is a very interesting 
case since, there will be an expansion fan at the change of cross section. This 
is an isentropic phenomenon [19]. As a result, there will not any entropy 
generation along the Wall. So, any entropy generation in the computations 
for this case will be purely numerical in nature. The Mach Contours obtained 
for this case are presented in fig.4.11. The Wall entropy along the Wall are 
shown in fig.4.12a. Though the Mach Contours donot show any difference 
between the two types of Wall boundary condition, the entropy plots highlights 
the difference. They clearly show that the entropy generation due to the 
Mirror Wall boundary condition is more as compared to that due to the Mixed 
Wall boundary condition, thus confirming the superiority of the Mixed Wall 
boundary condition over the Mirror Wall boundary conditions. 

It has been established that the strong formulation is less dissipiative [I], 
As stated in the earlier sections, the Mixed Wall boundary condition makes 
use of the strong formulation to determine the fictitious state. The strong 
formulation takes into account the curvature effect. Then the interior state 
is determined in a way similar to the Mirror condition. The mirror boundary 
condition, as we stated earlier, does not take into account the curvature effects 
and also it is dissipiative. The first order computations of Mirror conditions 
in case of finite volume, will have variables at the walls same as the cell val- 
ues themselves. Therefore, these wall boundary condition assumes the zero 
normal pressure and density gradients on the wall and thus violates Crocco’s 
relation and momentum equation. When a higher order reconstruction pro- 
cedure is adopted, though the spatial variations of the variables are properly 
represented, numerical experiments indicate that this alone is inadequate to 
make the Mirror boundary condition less dissipiative. 
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(a) Grid 1: Airfoil NACA0012 



(b) Grid2:Expansion channel 

Figure 4.1: Grids used 
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(a) Mirror Wall Boundary Condition 



(b) Mixed Wall Boundary Condition 


Figure 4.2: MACH CONTOURS WITH VAN LEER SCHEME ( Mach No. = 
0.63 and Alpha = 2.00 degrees) 
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(a) Mirror Wall Boundary Condition 



(b) Mixed Wall Boundary Condition 


Figure 4.3: WALL ENTROPY FOR THE VAN LEER SCHEME ( Mach No 
= 0.63 and Alpha = 2.00 degrees) 
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(b) Wall entropy with the van Leer scheme 

Figure 4.6: MIXTURE OF MIRROR AND MIXED BOUNDARY CONDI- 
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(a) Mirror Wall Boundary Condition 



(b) Mixed Wall Boundary Condition 

Figure 4.7: MACH CONTOURS WITH LEER SCHEME ( Mach No. = 0.8 
and Alpha =1.25 degrees) 
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(a) Mirror Wall Boundary Condition 



(b) Mixed Wall Boundary Condition 


Figure 4.8: WALL ENTROPY FOR THE LEER SCHEME ( Mach No = 0.8 
and Alpha = 1.25 degrees) 
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Figure 4.12: WALL ENTROPY COMPARISON 
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